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We have performed molecular dynamics simulations to study the efltect of an external electric field 
on a macroion in the solution of multivalent Z : 1 salt. To obtain plausible hydrodynamics of the 
medium, we explicitly make the simulation of many neutral particles along with ions. In a weak 
electric field, the macroion drifts together with the strongly adsorbed multivalent counterions along 
the electric field, in the direction proving inversion of the charge sign. The reversed mobility of the 
macroion is insensitive to the external field, and increases with salt ionic strength. The reversed 
mobility takes a maximal value at intermediate counterion valence. The motion of the macroion 
complex does not induce any fiow of the neutral solvent away from the macroion, which reveals 
screening of hydrodynamic interactions at short distances in electrolyte solutions. A very large 
electric field, comparable to the macroion unscreened field, disrupts charge inversion by stripping 
the adsorbed counterions off the macroion. 

PACS numbers; 61.25.Hq, 82.45.-h, 82.20.Wt 



I. INTRODUCTION 

The concept of electrostatic screening has been well 
known since the work by Debye and Hiickel of early 
20t/i century [Q. In recent years, screening by strongly 
charged ions was found to result in counterintuitive 
phenomena such as attraction between like-charged 
macroions m- |3|. |4|, 5(1, and inversion of macroion charges 

ICP'rl 

122, m 



i2|, p, |T1, @ p H, 

Charge inversion was studied by 
experiments using both colloids and biological materi- 
als IliR, Jl IJl^JllJll'Jil 111' analytical theories 
@,^|r£l,^3~iirM"^i^ and by Monte Carlo and 
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direct method to observe 
This 



molecular dynamics (MD) simulations |15 

Experimentally, the most direct method 
charged colloids and macroions is electrophoresis, 
method is also the prime candidate for the technique of 
observing charge inversion. Although straightforward, 
this approach involves many questions upon a closer look. 
For example, does the macroion drift along with its ad- 
sorbed multivalent ions when an external electric field is 
applied? How many multivalent ions are attached to the 
macroion strongly enough to drift together? How is the 
drift affected by the solvent viscosity and the counterflow 
of monovalent ions? What happens when the field be- 
comes very strong? What is the field strength sufficient 
to disrupt the charged complex? These are the funda- 
mental questions necessary to address in order to bridge 
theoretical concepts of charge inversion and experimental 
observations. 

It should be born in mind that electrophoresis in gen- 
eral has quite a few delicate aspects. Some peculiar ones 
attracted significant attentions recently j2^ (see also the 
review article and references therein). The electric 
field acts not only on the macroions, but on every ion in 
a solution. In many cases, this leads to effective screen- 



ing of hydrodynamic interactions which otherwise may 
be very significant. In the simulations reported below, 
we have actually observed such short-range screening of 
hydrodynamic interactions in the system comprising of a 
macroion, counterions and colons (Scc.IIIBl). 

One of the difficulties in simulating charge inversion 
under electrophoresis consists in subtle interactions of a 
macroion with surrounding ions and neutral solvent. A 
naive use of the Langevin equation, assuming that ev- 
ery ion (radius R) in the system is subject to Stokesian 
friction —Qi:rjRv and white noise random forces that bal- 
ance the friction through the fluctuation-dissipation the- 
orem, is not justifiable. A simple counterexample would 
be two closely located spheres. Since other particles (ei- 
ther ions or neutral solvent) are excluded from the vol- 
ume between spheres, neither their corresponding friction 
forces nor random forces add to each other. One way of 
going around this problem is to incorporate macroion- 
solvent interactions using the Oseen tensor psf . This is, 
however, not easy to implement in numerical simulations, 
because the interactions produce complicated spatial cor- 
relations among random forces. Therefore, we address 
this problem by a direct approach introducing explicit 
neutral particles to deal with the macroion-salt-solvcnt 
interactions in the molecular dynamics simulations. 

The explicit simulation of the solvent molecules is of 
course costly. In this paper, with the limited computa- 
tional resources, we restrict ourselves to the system with 
only the Z : 1 salt, and no 1 : 1 salt. It is needless to say 
that in real water solvent there is always some amount of 
1 : 1 salt. In this sense, our present paper demonstrates 
the principle that charge inversion is a phenomenon ob- 
servable by a direct electrophoresis experiment. Further 
study including the 1 : 1 salt will be required to compare 
results with realistic systems. Here, we will specify the 
deviations arising from the lack of the 1 : 1 salt. Also, we 
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will confine ourselves to the study of a single macroion 
interacting with surrounding salt ions and solvent. 

Our plan in this study is to examine electrophoresis of 
a spherically shaped, uniformly charged macroion. We 
will systematically measure the mobility of the drifting 
macroion complex placed under an external electric field 
by molecular dynamics simulations. We first show that 
charge inversion does take place in a solution containing 
multivalent counterions, as manifested by the inverse mo- 
bility under the weak electric field. We then look at the 
dependences of the reversed mobility on the parameters, 
such as the concentration of co- and counterions and the 
surface charge density of the macroion. Finally, we con- 
sider the strong electric field regime and show that the 
strong field disrupts the charge-inverted macroion com- 
plex and terminates the charge inversion phenomenon. 



II. SIMULATION MODEL AND PARAMETERS 
A. Description of the Model 

We adopt the following model, with a, e, and m be- 
ing the units of length, charge and mass, respectively. 
(We have in mind a ^ 2A and m ~ 40 a.m.u.) A 
macroion with negative charge Qo between — 15e and 
— 180e is surrounded by counterions of a positive 
charge Ze and N~ w 300 colons of a negative charge 
— e. The system is maintained in overall charge neutral- 
ity, Qo + N^Ze — N~e = 0, which determines N'^ for a 
given Z. The mass of the macroion is M = 200m, and 
the mass of the co- and counter-ions is m. We also in- 
clude TV* neutral particles with mass m/2, where we note 
the mass of water molecule against that of K+ or Ca^+ 
ions. Approximately one neutral particle is located in 
every volume element (1.5a)'^ « (SA)'^ inside the simula- 
tion domain, excluding the locations already occupied by 
the macroion and other ions, which typically yields 8000 
neutral particles. These particles are confined in a rect- 
angular box of size L, with periodic boundary conditions 
in all three directions. Most of the runs are performed 
for L — 32a, except one series of the runs intended to 
test the finite size effect of the domain (Fig.0) described 
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in Sec. Ill Bl 



In addition to the Coulomb forces, all particles inter- 
act through the repulsive Lennard- Jones potential (pLj ^ 
4£[(a/r,,)i2 - (a/ry)6] for = |r, - r,I < 2i/V, and 
(/>ij — — e otherwise. Here is the position vector of 
the i-th particle, and a is the sum of the radii of two 
interacting particles, which are chosen as follows: ra- 
dius of the macroion, Rq, is between 3a and 5a, coun- 
terions and colons have radius a, and neutral particles 
a/2. We relate e with the temperature by e = fcsT, 
and choose fc^T = e^/5ea (we assume spatially ho- 
mogeneous dielectric constant e). The Bjerrum length 
is thus Xb = e^/ekBT = 5a. For the parameters of 
this temperature, the valence Z = 3, and the num- 
ber of colons N~ = 300, the ionic strength becomes 



After knowing the ionic strength, one is tempted to 
compute the Debye length which turns out to be ^ 
0.45a. We should emphasize that this number does not 
have much meaning for the system under study, because 
we work in the domain very far from applicability of the 
standard Debye-Hiickel theory. In particular, the average 
number of ions in the volume X'jj turns out smaller than 
unity. This is by no means surprising, because there is 
no charge inversion in the Debye-Hiickel theory, and to 
examine charge inversion we need to go to the regime 
where this theory fails. 

Since we do not have the 1 : 1 salt in our simulation, we 
should keep in mind that correlations between strongly 
charged Z-ions may be significant even away from the 
macroion. Indeed, in the theory of charge inversion [ p^ , 
the role of correlations is emphasized for the Z-ions in 
the vicinity of the macroion, where their concentration 
is particularly large. In our system, the concentration of 
the Z-ions is not very small even in the bulk, and we face 
the situations in which correlations between the Z-ions 
away from the macroion affect our results. For such cases, 
we make additional runs with the reduced concentration 
in the Z : 1 salt. However, we do not use this reduced 
concentration for all the runs, because such system is 
more prone to noises and fluctuations, requiring larger 
statistics to obtain reliable results. 

Calculation of the Coulomb forces under the periodic 
boundary conditions involves the charge sum in the first 
Brillouin zone and their infinite mirror images (the Ewald 
sum[^). The sum is calculated with the use of the 
PPPM algorithm g We use (32)^ spatial meshes 
for the calculation of the reciprocal space contributions to 
the Coulomb forces, with the Ewald parameter a « 0.262 
and the real-space cutoff Vcut — Ri + 10a, where Ri is the 
radius of the i-th ion. A uniform electric field E is applied 
in the x-direction. 

When starting the molecular dynamics simulation run, 
we prepare an initial state by randomly positioning all 
the ions and neutral particles in the simulation domain 
and giving Maxwell-distributed random velocities corre- 
sponding to the temperature Tinitiai- We integrate the 
Newton equations of motion with the use of the leapfrog 
method |Q, which is equivalent to the Verlet algorithm. 
In the absence of the electric field (£■ = 0), our system 
is closed, and its energy is conserved. After an initial 
transient phase, the distribution of velocities relaxes to 
a Maxwellian, corresponding to an equilibrium sampling 
of the microcanonical ensemble. This new Maxwell dis- 
tribution has the temperature T, which is a little higher 
than Tinitiai, because of the release of the potential en- 
ergy due to screening, i.e., local balancing of charges. We 
adjust Tinitiai such that fc^T — e. This makes e to be the 
unique relevant scale of energy, and, accordingly, we put 
T = ay/m/e as the unit of time. We choose At — O.OIt 
as the integration time step. The simulation runs are 
executed up to lOOOr. 
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B. Hydrodynamic interactions, their screening, 
and the temperature control 

When an external electric field is present, it performs 
work on the system. The corresponding energy, which 
is the Joule heat, is transferred to background neutral 
particles through collisions with accelerated ions. In our 
work, we simulate an electrophoretic bath that is kept at 
a constant temperature T. For this purpose, we pretend 
that all neutral particles go through the thermal bath of 
infinite heat capacitance, whenever they cross the bound- 
aries of the simulation domain (at the center of which the 
macroion is located at every moment). Operationally, we 
refresh the velocities of the neutral particles according 
to the thermal distribution when they cross the domain 
boundaries. This procedure maintains temperature sta- 
bly to within 5%. 

Two closely related factors are potentially dangerous 
as they might affect the molecular dynamics simulation 
results. One factor is the finite-size of the simulation 
domain, and the other is the long-range character of 
both hydrodynamic and Coulomb interactions. These 
problems become particularly important because some of 
the methods to simulate a constant temperature thermal 
bath are believed to lead to the screening of hydrody- 
namic interactions. This is obviously not acceptable in 
the simulation where the long range character of hydro- 
dynamic interactions is expected to be important | |3^ . 

Following (26[ ^ , we argue that hydrodynamic inter- 
actions are in fact effectively screened in our system and, 
therefore, that the domain size in our simulation is quite 
sufficient and the heat bath procedure described above is 
benign and reliable. 

To understand the situation, it is worth discussing the 
major point - the screening of hydrodynamic interactions. 
To begin with, why are hydrodynamic interactions long 
ranged? That fact can be understood well from the point 
of view of momentum conservation. Consider a particle 
immersed in a fluid and suppose that we pull this par- 
ticle with an externally applied force (such as gravity). 
Obviously, this force transfers momentum into the sys- 
tem and, however large the container may be, this mo- 
mentum must be transported away through the container 
walls. This necessitates the long range character of the 
hydrodynamic field. More accurately, if we surround our 
object by an arbitrarily large closed surface, then (un- 
der the stationary conditions) the outflow of momentum 
through this surface must be equal to the inflow of mo- 
mentum due to the external force. Because of the obvi- 
ous analogy with the Gauss theorem in electrostatics, we 
see that hydrodynamic field is just as long ranged as the 
Coulomb field (even though it has a more complicated 
vector structure). 

The above description must be modified significantly 
when the applied external force is due to the electric field 
and the solution is neutral as a whole. In this case, there 
is no overall inflow of momentum into the system, and 
therefore, there should not be any outflow through the 



walls of the container. More specifically, if there is one 
negatively charged macroion as in our simulation, it is 
surrounded by positively charged counterions and nega- 
tively charged colons such that the total charge of the 
crowd effectively vanishes at some finite distance. In the 
simple case of the Debye-Hiickel theory, this happens at 
about the Debye screening length A d from the macroion 
surface. Therefore, no momentum is transported further 
away, and hydrodynamic interactions are screened at the 
distances of the order of ^ . 

In this paper, we treat a more complicated situation 
in which the Debye-Hiickel theory does not apply and 
it is not easy to judge a priori at which distances the 
hydrodynamic interactions are screened. Nevertheless, 
since the system is neutral, hydrodynamic interactions 
must be scree ned. We therefore perform special test (de- 
scribed in Sec. Ill B 1 ) looking at the dependence of the 
macroion drift on the simulation domain size. We find 
that the drift is essentially size-independent at L > 20a 
which is the direct manifestation of the screening of hy- 
drodynamic interactions. 

Since hydrodynamic interactions are screened, our sim- 
ulation is not very sensitive to the method of maintaining 
the constant temperature. To test it, we have performed 
separate runs at weak electric fields, E < O.Se/ae, in 
which case we can run the simulation even without any 
heat drainage for a significant period of time before any 
noticeable heating of the system; the results of these con- 
trol runs are within error bars of the data obtained using 
the thermal bath (Fig||). 



III. SIMULATION RESULTS 
A. General Properties 

Our simulation results are shown in Figs. 1-7. Fig. |^ is 
a bird's-eye view of (a) all the ions and (b) the vicinity of 
the macroion. Counterions are shown in light blue and 
colons in dark blue (neutral atoms are not shown) . In this 
figure, the macroion charge is taken to be Qq = — 30e, 
its radius Rq — 3a, counterion valence Z ~ 3, and the 
electric field E = 0.3e/ea. It is seen that the macroion 
is predominantly covered by the counterions. As in the 
case without the electric field , the radially integrated 
charge has a sharp positive peak at a distance about a 
from the macroion surface. This peak is due to the posi- 
tive counterions being adsorbed on the macroion surface. 
The value of the peak charge under the conditions of 
Fig.| is Qpeak ~ 1.6|Qo|. 

Fig. 1^ demonstrates the time history of (a) the "peak" 
charge and (b) the macroion drift speed for the parame- 
ters of Figj^. At time t = lOr, we switch on the external 
electric field. There is a short transient phase during 
which a charge-inverted complex is formed through ad- 
sorption of counterions to the macroion and condensation 
of colons on the counterions. This process is reflected in a 
rather quick rise in Qpcak, as is shown in Fig.H(a). After 
the transient phase, we observe a drift of the macroion 
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FIG. 1: Bird's-eye view of (a) all the ions in the simulation 
domain and (b) the screening ion atmosphere within 3a from 
the macroion surface. A macroion with charge Qo = — 30e 
and radius Ro = 3a is a large sphere in the middle; counte- 
rions {Z = 3) and monovalent colons are shown by light and 
dark blue spheres, respectively. The arrow to the right shows 
the direction of the electric field (x-axis), with E = O.Se/ea. 

in the positive direction along the applied field. The fact 
that the drift velocity is positive for the negative bare 
charge of the macroion (Qo < 0) is a direct manifestation 
of charge inversion such that coimterions are so strongly 
bound that they pull the macroion with their motion. 

Note that the drift velocity shown in Fig.||(b) is small 
compared to the thermal velocity vq of neutral particles, 
(Vx) ~ O.OSwq. Under this condition, exchange of mo- 
mentum between the macroion and neutral particles is 
slow, and it requires many collisions (compare the sim- 
ilar system in Ref . ||3^ ) . Therefore, in terms of hydro- 
dynamics, we are in the linear regime. It means that 
friction force should be linear in the macroion velocity 
and we expect the average drift speed to be given by the 
force balance condition, Q*E — vVx — 0, where Q* is 
the effective net charge of the macroion complex and v 
is the hydrodynamic friction coefficient. We shall discuss 
later the possibilities of determination of the effective net 
charge Q* based on this condition. 

Fig. ^ also shows significant temporal fluctuations in 
the drift speed. Inspection reveals that they are larger 
than what one expects for random kicks of neutral parti- 
cles. These fluctuations indicate that neither the counte- 
rions permanently stick to fixed points on the macroion 
nor the colons attach to the counterions, but that they 
are being replaced from time to time. The fluctuations 
of this type are actually seen in a video movie. 

B. Parameter Dependences 

1. Linear Regime 

The dependence of the average macroion drift speed 
Vdrift on the electric field is shown in Fig.||. In this fig- 
ure, we show together the results of several runs, corre- 
sponding to different values of the macroion charge Qo 
and macroion radius i?o. First and foremost, the sign of 
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FIG. 2: Time history of (a) the peak charge Qpeak (de- 
fined as the maximum of radial charge distribution around 
the macroion center) and (b) the macroion speed Vx normal- 
ized by thermal velocity of neutral particles wo . The macroion 
complex drifts positively along the external electric field of 
-E > 0, which directly indicates the inversion of the charge 
sign. 



the drift velocity in moderate fields corresponds to the 
sign of inverted charge. This is the central observation 
of our work. The figure clearly demonstrates the overall 
pattern of the drift velocity dependence on the applied 
field, beginning with the linear regime in a weak field, 
followed eventually by a breakdown of charge inversion 
in a strong field. 

Let us discuss the linear drift regime for small elec- 
tric fields E < 0.2|Qo|/e-Ro' where Vdrift increases lin- 
early with the field strength. This regime corresponds to 
the usual Ohm's law, where the net charge of the com- 
plex is insensitive to the strength of the electric field. A 
macroion drifts together with its strongly adsorbed coun- 
terions as a complex. That is, the electric field is not 
strong enough to affect the binding of counterions to the 
macroion. 

At this stage, it is necessary to check the issue of hy- 
drodynamic interactions and their screening. For that 
purpose, we show in Fig. ^ the effect of the simulation 
domain size L on the macroion drift speed. By a series 
of the runs with different domain sizes and under fixed 
number density of neutral particles and ionic strength, 
we have confirmed that ed L — 32a, which is the domain 
size for the majority of our simulations, the domain-size 
dependence is essentially leveled off. This is to be con- 
trasted with polymer chains p3| , in which case hydrody- 
namic interactions lead to the finite size effect essentially 
linear in a/L. 

We have also inspected the fluid flow of neutral parti- 
cles around the macroion. When rapid fluctuations are 
averaged out, this flow field does not exhibit any patterns 
protruding away from the macroion. This fact implies 
that the fiow of neutral particles induced by the motion 
of the macroion and other ions is screened at short dis- 
tances |2^, . 
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FIG. 3: Dependence of the macroion drift speed Vdrift (in 
the units of vo , the thermal speed of neutral particles) on the 
electric field E for a macroion of various radii and charges: 
Ro = 3a and Qo = — 30e (filled circles); i?o = 4a and Qo = 
— 50e (open triangles); and 7?o = 5a and Qo = — 80e (open 
circles); Rq — 5a and Qo = — 51e (open squares). The valence 
of counterions is Z — 3. 
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FIG. 4: Effect of the finite domain size L on the drift speed 
Vdrift for a macroion with Ro — 3a, Qo = — 30e, the electric 
field E = 0.3e/ae, and counterion valence Z — 3. The error 
bars correspond to the root mean square velocity fluctuations 
as seen in Fig. ^. 



The saturation of the a/L dependence (Fig. Q) and the 
inspection of neutral particle flow patterns both confirm 
that hydrodynamic interactions are screened in our sys- 
tem, thus making reliable our simulation approach based 
on the finite domain and the heat bath. The thermal 
bath at the distant boundaries does not affect the mea- 
sured macroion mobility. 

The small electric field regime is characterized by the 
mobility, ^ — (Vx) /E. This quantity is plotted in Fig. ^ 
as a function of the macroion bare charge Qq, or, more 
specifically, on the surface charge density of the macroion 



2. Can we measure net charge based on the mobility 
measurement data? 



Let us now discuss a practically important question: 
Can we determine the effective net charge of the macroion 
complex Q* based on the data of mobility measurements, 
i.e., based on the data of Fig I - Fig.'l ? Physically, 
as we have already mentioned, the charge Q* is deter- 
mined by the force balance condition Q*E — vVx — 0, 
or Q* = ^v. Therefore, determination of the net charge 
requires knowledge of both the friction coefficient v and 
the mobility fi. 

Importantly, the friction coefficient cannot be deter- 
mined by the usual Stokes formula I's — STrryi?, where 
?7 is the solvent viscosity. The problem is that the real 
friction is enhanced by the screening of hydrodynamic 
interactions [^6[ The Stokes formula is simply un- 
derstood by that the friction force in general should be 
proportional to ri{v/A)B'^, where A is the length scale 
over which velocity changes, is the relevant surface 
area. For the Stokes problem, we have A ^ B ^ Rq. By 
contrast, when the drift is caused by the action of the 
electric field on the overall neutral system, the distance 
A becomes much smaller. In the Debye-Hiickel theory, it 
turns out to be of the order of the Debye screening length, 
A Xd- In this case, the friction coefficient becomes 



V = vsRo/\d 



(1) 



i.e., it is enhanced by the factor Rq/Xd compared to the 
usual Stokesian friction. Although we work under the 
conditions where the Debye-Hiickel theory is not appli- 
cable, and we do not know exactly which length should 
be there instead of Ad in Eq.(P, this length is clearly 
smaller than Rq and independent of Rq. Therefore, the 
friction coefficient is proportional to Rq - unlike more 
familiar Stokes case where it scales like Rq. 

Under usual circumstances where the charge Q* is 
known, the friction coefficient scaling as Rq implies that 
mobility fi = Q* /v cc Q* / Rq is determined by the sur- 
face charge density, not by the charge and the surface 
area separately. This fact was known to M.Smoluchowski 
already a century ago [^5| . In our simulation, effective 
charge is not known a priori, and the logic needs to be re- 
versed. Figure |^ indicates that mobility /j, is essentially a 
constant when the macroion bare surface charge density 
is not too small (|(5o|a^/ei?o > 3). Given that ly (x Rq, 
we conclude that the effective bare charge Q* — is 
proportional to the surface area of the macroion; namely, 
charge inversion is characterized by the overcharging den- 
sity. This agrees with the theory ^ . 

In order to perform at least very rough quantitative es- 
timate of charge Q* based on the mobility data, we need 
to know the Stokesian friction coefficient vs (or equiva- 
lently, we need to know the viscosity of our model sol- 
vent). We measure it in a separate molecular dynamics 
run, by observing an exponential decay of the macroion 
velocity starting from 0.5wo for the case without an elec- 
tric field. We find i's ~ 9.3m/r for a spherical particle of 
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the radius Rq = 3a and vs ~ 18.2m /t for the macroion 
complex with adsorbed counterions and coions. These 
two estimates provide lower and upper bounds for the 
charge Q*. Assuming Rq/Xo « 6 and « 0.5/io (sat- 
uration regime in Figj^) yields the inverted charge Q* 
between 7e and 20e. This is in rough agreement with the 
Qpeak measurements. 

Special attention must be paid to the small bare 
surface charge density case, for which Fig.^] indicates 
the decrease in the reversed mobility. For some cases, 
the reversed mobility even disappears altogether, chang- 
ing to normal, non-reversed mobility, /i < when the 
macroion bare surface charge density decreases to about 
IQala^/eRg ^ 1. It turns out that this is the man- 
ifestation of correlations between Z-ions in the bulk 
solution away from the macroion. Indeed, when the 
macroion is only weakly charged, the correlations of Z- 
ions in its vicinity are not much stronger than in the 
bulk, which suppresses charge inversion. Simple estimate 
shows that, under the conditions when fi gets small or 
negative in Fig.||, the "Wigner cell" radius of the Z-ions 
on the macroion surface, Rw — 2i?o(eZ/|Qo|)^^^ (which 
is about 3.4a for IQola"^ /cRq ~ 1 and Z = 3) becomes 
comparable to the average spacing between the .Z-ions in 
the bulk. The inspection of the radial charge distribu- 
tion functions around the macroion for this case agrees 
with this interpretation. It shows that the counterions 
are loosely bound to the macroion, and that the coions 
form relatively strong bonds with the counterions and 
drift together with them. 

To examine the above interpretation further, we per- 
form special runs with reduced concentration of the Z : 1 
salt. The results of these runs are shown in Fig.|| with 
filled circles and triangles for the iV+ = (A^^e-I- \ Qo\)/Ze 
Z-ions with N~ = 90 and 30 negative coions, re- 
spectively. As anticipated, with fewer Z-ions in the 
bulk, charge inversion is not interrupted even at small 
macroion bare charges. We regard this a convincing proof 
that charge inversion occurs at small concentration of the 
Z-ions despite a larger entropy penalty. 

Figure || also shows with crosses the data of the con- 
trol run performed under the condition of the weak elec- 
tric field without any heat drainage (see Sec. II B). These 
data are within error bars of the cases with the simulated 
thermal bath. 

The dependence of the macroion mobility fj, on the 
valence of the counterions Z in Fig|6| is physically inter- 
esting, and also important for application purposes. For 
the cases shown with filled and open circles, the surface 
charge density of the macroion is chosen nearly the same 
|Qo|/-Rn ~ 3 so that they reside in the saturation regime 
of Fig.H. We emphasize that the mobility for these cases 
is given by the same curve. The mobility of the macroion 
is negative when counterions are monovalent Z = 1, be- 
cause there is no charge inversion but only regular Debye 
screening. Thus, the charge inversion phenomenon does 
not occur in the solution of monovalent salt (provided 
that the co- and counter- ions have the equal radius). For 




FIG. 5: Dependence of the macroion mobility /i on the sur- 
face charge density Qo/Ro for the macroion radius _Ro = 3a 
(open triangles), Ro = 5a (open circles), and Ro — 7a (open 
diamonds), where fio ~ ^^o/dQo"' |/e(-Ro''^)^) with q'^'' = 
— 30e and R'^^ = 3a. The valence of the counterions is Z — S, 
the number of the Z : 1 salt is iV+ = {N~e + \Qo\)/Ze and 
N~ — 300, the electric field is i5 = 0.3e/ae, and the temper- 
ature is e^/eaksT = 5. The filled circles and triangles show 
the cases with reduced number of the Z : 1 salt such that 
= 90 and 30, respectively. The crosses are the reference 
data obtained without the thermal bath for Ro — 5a and the 
Z : 1 salt with iV" = 300. 



Z > 2, charge inversion does take place, as manifested 
by the positive mobility. These observations agree with 
a previous study for planar charged surfaces jl^. A re- 
markable finding is that the mobility here is maximized 
for the intermediate valence, Z 4, unlike the peak in- 
verted charge that accounts for static charge distribution 
of mainly counterions p3| . 

It is also noted in Fig]^ that, when the surface charge 
density is reduced, both the magnitude of reversed (pos- 
itive) mobility and the range of Z where it occurs shrink 
as shown by square symbols in the figure. The mobility 
for divalent Z-ions is now negative. This corresponds to 
the small surface charge density regime \Qo\/Rq ^ 1.6 in 
Fig.^. Yet, the mobility is maximized for the intermedi- 
ate valence, Z = 5 in this case. Also shown in FigJ^ are 
the results of the control runs performed with reduced 
number of Z-ions (discussed above in connection with 
Fig. H). They again reproduce the optimal valences for 
charge inversion. 

Although somewhat speculative, we can try to apply 
the result of Fig.^ to explain the electrophoretic mobil- 
ity measurements of nucleosome core particles in cation 
solution iQ. What was observed is the increase in the 
magnitude and range of the cation concentration for oc- 
currence of reversed mobility when spermidine salt (-1-3) 
was replaced with spermine salt (+4), while charge rever- 
sal was not observed for any concentration of Mg"*"^. It 
is worth stressing that nucleosome particles is a compli- 
cated system in which many aspects may be important. 
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FIG. 6: Dependence of the macroion mobility /i on the 
valence of the counterions Z for the runs: Ro = 3a and 
Qo ~ — 30e (filled circles), Ro — 5a and Qo = — 80e (open 
circles), and Ro = 5a and Qo ~ — 40e (filled squares). Here, 
the number of the Z : 1 salt is iV+ = {N~e + \Qo\)/Ze and 
A'^" « 300, the external electric field is i? = 0.3e/ae, where 



MO = wo/(|Qj,"^|/e(-R^°')^ with Ql,"' = -30e and R'^' = 3a. 
A series of the runs with reduced number of the Z : 1 salt 
~ 30, Rq = 5a, and Qo ~ — 80e are shown with open 
triangles. 
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FIG. 7: Dependence of the macroion mobility fi on the salt 
ionic strength, nj — {Z^N^ + N^^jlJ" . The parameters are 
Qo = — 30e, Ro — 3a, E — 0.3e/ae, Z — 3 and e^/eaksT — 5, 
which yield the Debye length Ad ~ 0.45a for nj = 0.04a~"^. 



reversed mobility at small ionic strength in Fig.^ agrees 
with the colloidal mobility measurement for the case with 
trivalent salt LaCls [|[ |6|. 



What we would like to say here is that our results may 
be at least one of the factors relevant to the experiments 
reported in the literature Q. 

The dependence of the macroion mobility on the salt 
ionic strength, m = {Z'^N^ + N~)/L^ is shown in Fig.|^ 
for the counterion valence Z — 3 and the temperature 

/taksT — 5. Here, the ionic strength is related to the 
Debye length by Ad = {eksT /d.Tmie'^Y/'^ . The mobility 
increases quite rapidly for small ionic strength, and is 

1/2 

well fit by /Lt oc as shown by a dashed curve. 

It is a legitimate concern to ask whether the data of 
Fig.0 are affected by the correlations of Z-ions in the bulk 
which we discussed in connection with Fig. ^ The an- 
swer is negative jthese data correspond to the saturation 
regime of Fig. ^. The distance between Z-ions in the 
bulk drops to the value comparable to the Wigner-Seitz 
cell only when nia^ gets as large as 0.1 or more. 

The increase in the mobility with ionic strength con- 
tradicts the intuition based on the Debye-Hiickel theory, 
and deserves a comment. As it is explained in detail in 
psf and understood by a number of authors cited there, 
charge inversion itself grows with ionic strength. This 
happens because charge inversion is the result of inter- 
play between the repulsion of counterions from each other 
and the attraction of each of them to its own correlation 
hole. The latter occurs at a much shorter distance than 
the former, and only the repulsion is strongly affected by 
screening. This is why the amount of charge inversion, 
hence the macroion mobility, grows with increasing ionic 
strength. 

It is also worth mentioning that the quick rise in the 



3. Nonlinear Regime 

Let us now return to Figj| to discuss the regime that 
is nonlinear in the applied electric field. As the figure 
indicates, the charge-inverted shell around the macroion 
is destroyed for large electric fields. Moreover, the criti- 
cal field Ec at which this happens is independent of the 
macroion size, which leads us to an estimate 

Ec « 0.5|Qo|/ei?g . (2) 

This result is quite interesting. Indeed, \Qf)\/eR^ is the 
electric field on the macroion surface produced by the 
macroion bare charge. Why does the critical field scale 
with the bare charge of the macroion instead of the net 
charge of the complex ? The reason is due to correla- 
tions between screening ions. We noted while discussing 
Fig.^ that the counterions on the macroion surface are 
being replaced from time to time. Consider how one Z- 
ion can depart from the macroion surface. Since this 
ion is surrounded by a correlation hole on the surface, 
its departure requires work against the unscreened bare 
electric field of the macroion as long as its distance from 
the surface is smaller than the distance between the ad- 
sorbed Z-ions. Therefore, departure from the surface 
becomes possible when the external field becomes com- 
parable with this unscreened field; the charge-inverted 
complex is no longer stable at such a field strength. 

The critical electric field in realistic situations is esti- 
mated to be very large. For the parameters i?o ~ 20Aand 
Qo ~ 30e, the critical electric field becomes as large as 
Ec ~ 0.5(3o/e-Ro ~ QlV/fim, where we take into account 
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the dielectric constant of water e « 80 |36 . Although 
the critical field is large, it gives small energy to the elec- 
tric dipole of a water molecule, d « 2 x lO^^^esu ■ cm: 
Ecd/kBT ~ 0.11 < 1. This verifies the use of the model 
solvent of neutral particles in the present molecular dy- 
namics simulations. In practice, the applied electric field 
is not expected to disrupt the charge-inverted macroion 
complex. 



IV. SUMMARY 

In this paper, we performed molecular dynamics simu- 
lations with the use of neutral-particle solvent, and mea- 
sured the drift speed of a macroion to obtain its mobility 
under electrophoresis in a multivalent Z : 1 salt solution. 

A weak electric field pulled the macroion complex in 
the direction determined by the net inverted charge, in- 
stead of disrupting it. The reversed mobility of the com- 
plex, /i — Vdrift/E, was shown to be nearly constant for 
the weak electric fields. We showed the functional depen- 
dences of mobility in Figs. 3, 5 and 6 of Sec. Ill, respec- 
tively, against the electric field strength E, the surface 
charge density of the macroion Qq/Rq, and the valence 
Z of the counterions. The mobility was a function of the 
surface charge density, instead of the bare charge and 
radius of the macroion separately. The reversed mobil- 
ity increased rapidly with the salt ionic strength nj as 

1/2 

/i cx n/ . Interestingly, the reversed mobility took a 
maximal value at the intermediate valence of the coun- 
terions Z = 4. 

We confirmed the screening of hydrodynamic inter- 



actions at a few Debye length. No specific flow pat- 
terns of neutral particles, which one would expect for the 
sphere moving in a viscous fluid, were detected around a 
macroion. 

In the large field regime, although academic because 
of its huge value, electrophoresis was strongly nonlin- 
ear, and the field stripped the screening counterions off 
the macroion. The mobility of the macroion complex 
dropped significantly from that of the linear regime, and 
the sign of the mobility flipped back to non-reversed one 
above the critical electric field, which was approximately 
half the macroion unscreened field. 

In this study, we explicitly simulated the neutral par- 
ticles of the solvent to produce a reliable hydrodynamic 
background. On the other hand, the limits of computa- 
tional resources prevented us from inclusion of the 1 : 1 
salt. The screening in our system was exclusively ac- 
complished by the Z : 1 salt. For this reason, we are 
not ready to make a quantitative comparison of our data 
with the real experiments. The simulation including the 
1 : 1 salt is currently under way. 
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